#include "first_macro_2d.h"  /* Where global parameters are defined */
#include "structs_2d.h"      /* Where structures are defined */
#include "derivs_2d.h"       /* Where derivatives functions are defined */
#include "gen_2d.h"

/***********************   Global variables   ****************************/


/* ---------------------------  rkc -------------------------------------*/

struct field_array dv_temp;        /* derivative temporal value         */
struct field_array dv_sum;         /* sum values of deriv               */
struct field_array v_temp;         /* intermediate value of v           */

struct field_array *dv_temp_ptr = &dv_temp;    /* derivative temporal value        */
struct field_array *dv_sum_ptr = &dv_sum;      /* sum          values of deriv     */
struct field_array *v_temp_ptr = &v_temp;      /* intermediate value of v          */



/* ---------------------------- wave -----------------------------------*/

/* temporal-space first derivative values */

struct deriv_array dfields;   


/* auxiliary field values */


struct aux_array auxfields;

/* stencil fields */

struct s_field_array sfield_recv;
struct s_field_array sfield_send;


/* MPI */

MPI_Comm global;

/********************************************************************/



int main(int argc, char *argv[]) {

  /* variable declarations */

/*   int argc; */
/*   char *argv[]; */

  struct GRID_PAR grd;
  struct GRID_PAR *grd_ptr = &grd;
  FLOAT h;
  struct field_array y; 
  FLOAT swap;
  #ifdef ENERGY
  FLOAT E,W;
  #endif
  


  /*   Ploting names */

  struct PLOT_PAR plot; 
  struct PLOT_PAR *plot_ptr = &plot; 


  /*   Initial data parameters */

  struct INI_PAR init_parameters;
  struct INI_PAR *init_par_ptr = &init_parameters;
  
  
  /*   Function parameters  */


  struct FUNCTION_PAR equation_parameters;
  struct FUNCTION_PAR *equation_par_ptr = &equation_parameters;


  /* ---------------------------------------- Start MPI ------------------------------------------------- */

#ifdef MPI

  MPI_Status status;

  MPI_Init(&argc, &argv);

  MPI_Comm_dup(MPI_COMM_WORLD, &global);


 MPI_Comm_size(global, &grd.n_grids);
 if (grd.n_grids != N_GRIDS)  {
                             printf("Number of processes and of grids do not coincide");
                            MPI_Finalize();
	       return(1);
 }

 MPI_Comm_rank(global, &grd.grid);

#else 
grd.grid=0; 
#endif

  /* Parameters coming from first_macro */


#ifdef FILE_INPUT
  FILE *file_data_ptr;
#endif

  /* Get data from web page or data file */
 
#ifdef DEBUG_INPUT
printf("going into input function\n");
#endif 


INPUT_FUNCTION(grd_ptr, equation_par_ptr,init_par_ptr, plot_ptr);

#ifdef DEBUG_INPUT
 printf("out of input function\n");
#endif 

file_data_ptr = plot.input_data_file_ptr;
PRINT_MACRO_VALUE(RKX)
PRINT_MACRO_VALUE(DERIV)
PRINT_MACRO_VALUE(FF)
PRINT_MACRO_VALUE(ADISCO)
#ifdef DISSIPATION
PRINT_MACRO_VALUE(DISS)
#endif

  /* Allocate memory */ 

//exit(0);

#include "m_alloc.c"

/* the grid number and the number of grids */



  /* -------------------------- Grid specific parts --------------------------------- */

  plot.n_plots = N_PLOTS;
  
/* Array of names for plots */

 sprintf(plot.name[0],"M_A_N");
 sprintf(plot.name[1],"M_A_T");
 sprintf(plot.name[2],"M_A_X");
 sprintf(plot.name[3],"M_A_Y");
 sprintf(plot.name[4],"M_B_N");
 sprintf(plot.name[5],"M_B_T");
 sprintf(plot.name[6],"M_B_X");
 sprintf(plot.name[7],"M_B_Y");
 sprintf(plot.name[8],"M_C_N");
 sprintf(plot.name[9],"M_C_T");
 sprintf(plot.name[10],"M_C_X");
 sprintf(plot.name[11],"M_C_Y");

// Create the data file names
	{int np;
 		for (np=0; np <  plot.n_plots; np++) {
  		sprintf(plot.window_name[np], "%s_%d_%s_%d_N_%dx%d",plot.output_file_name,plot.prev_time_stamp,plot.name[np],grd.grid,grd.n_gridpts_1,grd.n_gridpts_2);
		}
	}
	
  plot.initial_x= grd.initial_x;
  plot.final_x = grd.final_x;
  plot.initial_y = grd.initial_y;
  plot.final_y = grd.final_y;



#ifdef DEBUG
{int np;
 for (np=0; np <  plot.n_plots; np++) {
   printf("Plotting name is: %s \n",plot.window_name[np]);
 }
}
#endif




/* Relation between fields and plot names */

  plot.pointers[0] = M_A;  
  plot.pointers[1] = M_A_T;
  plot.pointers[2] = M_A_X;
  plot.pointers[3] = M_A_Y;  
  plot.pointers[4] = M_B;  
  plot.pointers[5] = M_B_T;
  plot.pointers[6] = M_B_X;
  plot.pointers[7] = M_B_Y;  
  plot.pointers[8] = M_C;  
  plot.pointers[9] = M_C_T;
  plot.pointers[10] = M_C_X;
  plot.pointers[11] = M_C_Y;  

/* Number of gridpoints in plots */

/*   plot.grid_plot_pts_1 = N_GRID_PLOT_PTS_1;   */
/*   plot.grid_plot_pts_2 = N_GRID_PLOT_PTS_2; */

    
  /* Factor between gridpts and gridplotpts */
 
/*   plot.factor_1 = FACTOR_1;  */
/*   plot.factor_2 = FACTOR_2;  */

  /* Initial/Final value for time coordinate */

  plot.initial_time = grd.initial_time;       
  plot.final_time = grd.final_time;   

  /* exiting */

//MPI_Comm_free(&global);
//MPI_Finalize();
//fflush(stdout);
//exit(0);  

/* Open output file (some times used only for compatibility) */
/* printf("opening file\n"); */

  plot_ptr = ADISCO('O',  &plot, &grd, &y); 

  /* makes initial time-interval */

  h = (grd.final_time - grd.initial_time)/(FLOAT)(grd.data_steps*grd.int_steps);

  /* sends input data to file / browser */


if(grd.grid ==0){

	fprintf(OUT_PUT,"<li> Total number of Time Steps = %f </br>\n",(double)(grd.data_steps*grd.int_steps));

	fprintf(OUT_PUT,"<li> Number of Time Steps per unit time = %f </br>\n",1.0/h);

	fprintf(OUT_PUT,"<li> Time_Step / Space_Step_x = (h/(xf-xi)*n_gridpts) = %f </br>\n",h*(double)grd.n_gridpts_1/(grd.final_x-grd.initial_x));
	fprintf(OUT_PUT,"<li> Time_Step / Space_Step_y = (h/(yf-yi)*n_gridpts) = %f </br>\n",h*(double)grd.n_gridpts_2/(grd.final_y-grd.initial_y));
	fprintf(OUT_PUT,"</ul>%c</br>",10);
	
	fflush(stdout);
	
	/* send input data to the screen */
	
	#ifdef FILE_INPUT
	printf("Total number of Time Steps = %f \n",(double)(grd.data_steps*grd.int_steps));
	
	printf("Number of Time Steps per unit time = %f \n",1.0/h);
	
	printf("Time_Step / Space_Step_x = (h/(xf-xi)*n_gridpts) = %f \n",h*(double)(grd.n_gridpts_1-1)/(grd.final_x-grd.initial_x));
	printf("Time_Step / Space_Step_y = (h/(yf-yi)*n_gridpts) = %f \n",h*(double)(grd.n_gridpts_2-1)/(grd.final_y-grd.initial_y));
	printf("\n");
	
	fflush(stdout);
	#endif
}

  /*     creates initial data                            */

inidat(&y,grd_ptr,&init_parameters);

/* write initial data to file                          */ 

/* plot data */
      
      
		if(plot.time_stamp==0){
/* plot data */
  			plot_ptr = ADISCO('P',  &plot, &grd, &y); 
		}
		else{
// get data from files and re-write it on initial data
  			plot_ptr = ADISCO('R',  &plot, &grd, &y); 

// rename the data files
				{int np;
 					for (np=0; np <  plot.n_plots; np++) {
  					sprintf(plot.window_name[np], "%s_%d_%s_%d_N_%dx%d",plot.output_file_name,plot.time_stamp,plot.name[np],grd.grid,grd.n_gridpts_1,grd.n_gridpts_2);
					}
				}
		}

fflush(stdout);





// compute de energy and send it to output

#ifdef ENERGY
E = Energy(grd_ptr,&y,dv_temp_ptr,equation_par_ptr);
W = Winding_number(grd_ptr,&y,dv_temp_ptr,equation_par_ptr);

	if(grd.grid ==0){
			printf("time= %f, E= %f, W= %f \n",y.time,E,W);
		fflush(stdout);
	}	
#endif



/* stop program after making a mask in inidat */

#ifdef MAKE_MASK 
 exit(0); 
#endif


//exit(0);  


/* inipot(pot_ptr,&pot_parameters, &grd, plot_ptr); */


 


  /* Take data_steps */

  {long int k_outer;
  for (k_outer=1; k_outer<= grd.data_steps; k_outer++) {

#ifdef DEBUG
      printf("h = %f\n",h); 
      printf("time = %f\n",y.time);
#endif

integ(&y,h,grd.int_steps,grd_ptr,equation_par_ptr,FF,RKX); 


/* 	printf("time after integ in main = %f",y.a.time); */
/* printf("Out of integ \n"); */ 



E = Energy(grd_ptr,&y,dv_temp_ptr,equation_par_ptr);
W= Winding_number(grd_ptr,&y,dv_temp_ptr,equation_par_ptr);


#ifdef MPI
	if(grd.grid ==0){
#else
{
#endif 
#ifdef ENERGY 
printf("time= %f, E= %f, W= %f \n",y.time,E,W); 
fflush(stdout);
#else
printf("%d",k_outer);
#endif

fflush(stdout);
}





      /* Save values */


      plot_ptr = ADISCO('P', &plot, &grd, &y);  

/* printf("�</br>\n"); */  
/* printf("%d/%d, ",k_outer, grd.data_steps);  */
/* printf("%c",7); */ 
fflush(stdout);




  }



  }

#ifdef MPI

	if(grd.grid ==0){
#else
{
#endif 


fprintf(OUT_PUT,"<ul>%c</br>",10);
fprintf(OUT_PUT,"<li> Execution time = %u secs. ", (unsigned)(clock()/CLOCKS_PER_SEC));
fprintf(OUT_PUT,"</ul>%c</br>",10);

#ifdef FILE_INPUT
printf("\n");
printf("n%d  Execution time = %u secs. ", grd.grid, (unsigned)(clock()/CLOCKS_PER_SEC));
printf("\n");
#endif
fflush(stdout);
}

/* close output file */
  plot_ptr = ADISCO('C',  &plot, &grd, &y); 

#ifdef FILE_INPUT
  fclose(plot.input_data_file_ptr);
#endif


/* Free memory */



{int n_fields;
 for (n_fields = 0; n_fields < N_FIELDS; n_fields++){
free(dv_temp.u[n_fields]);
free(dv_sum.u[n_fields]);
free(v_temp.u[n_fields]);
free(y.u[n_fields]);
 }
}
{int n_fields;
 for (n_fields = 0; n_fields < N_DERIVS; n_fields++){
free(dfields.du[n_fields]);
 }
}
{int n_fields;
 for (n_fields = 0; n_fields < N_AUX; n_fields++){
free(auxfields.u_aux[n_fields]);
 }
}   
{int field;
 for (field = 0; field < N_FIELDS_MPI; field++){
free(sfield_recv.s_x_0_u[field]);
free(sfield_recv.s_x_1_u[field]);
free(sfield_send.s_x_0_u[field]);
free(sfield_send.s_x_1_u[field]);
free(sfield_recv.s_y_0_u[field]);
free(sfield_recv.s_y_1_u[field]);
free(sfield_send.s_y_0_u[field]);
free(sfield_send.s_y_1_u[field]);
 }
}
{
free(plot.plot_field);
}


/* exiting */

MPI_Comm_free(&global);
MPI_Finalize();

/* printf("%c",7); */
printf("finishing \n");
return(0);
}

// ----------------------------------------------------------------------------------------

// --------------------- auxiliary functions ----------------------------------------------

// ----------------------------------------------------------------------------------------


FLOAT Energy(struct GRID_PAR *grid_ptr,
		struct field_array *fields_ptr, 
		struct field_array *derivs_ptr,
		struct FUNCTION_PAR *function_par) {


    int ni_1 = (*grid_ptr).start_grid_1; 
    int nf_1 = (*grid_ptr).final_grid_1;  
    int ni_2 = (*grid_ptr).start_grid_2;  
    int nf_2 = (*grid_ptr).final_grid_2; 
    int n_gridpts_1 = nf_1 - ni_1;
    int n_gridpts_2 = nf_2 - ni_2;
    int grid = (*grid_ptr).grid;
    int coorx = 1;    
    int coory = 2;


    FLOAT xi = (*grid_ptr).initial_x;
    FLOAT xf = (*grid_ptr).final_x;
    FLOAT yi = (*grid_ptr).initial_y;
    FLOAT yf = (*grid_ptr).final_y;

    FLOAT x, y, m_t_nn, nn_x, nn_y, nn_z;

    FLOAT one_dN1 = (xf-xi)/(FLOAT)(nf_1-ni_1-1);
    FLOAT one_dN2 = (yf-yi)/(FLOAT)(nf_2-ni_2-1);

    FLOAT sigma = (*function_par).sigma; /* For dissipation */

    FLOAT mm = (*function_par).mm;
	FLOAT E = 0.0;
	
	
	
    FLOAT g_up_xx, g_up_yy, g_up_xy, g_up_yx, detg_1_2;
	FLOAT m_a_up_x, m_b_up_x, m_c_up_x, m_a_up_y, m_b_up_y, m_c_up_y;
	FLOAT norm_m_2;
	
	FLOAT m_a=0.0;
    FLOAT m_a_t=0.0;
    FLOAT m_a_x=0.0;
    FLOAT m_a_y=0.0;

    FLOAT m_b=0.0;
    FLOAT m_b_t=0.0;
    FLOAT m_b_x=0.0;
    FLOAT m_b_y=0.0;

    FLOAT m_c=0.0;
    FLOAT m_c_t=0.0;
    FLOAT m_c_x=0.0;
    FLOAT m_c_y=0.0;

 {register int grid_ind1, grid_ind2;
  FLOAT Eint;
 
 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){
   for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){

         g_up_xx = auxfields.u_aux[G_UP_XX]GRD_IND;
         g_up_yy = auxfields.u_aux[G_UP_YY]GRD_IND;
         g_up_xy = auxfields.u_aux[G_UP_XY]GRD_IND;

	 detg_1_2 = 1./sqrt(g_up_xx*g_up_yy - g_up_xy*g_up_xy);


       m_a   = (*fields_ptr).u[M_A  ]GRD_IND;
       m_a_t = (*fields_ptr).u[M_A_T]GRD_IND;
       m_a_x = (*fields_ptr).u[M_A_X]GRD_IND;
       m_a_y = (*fields_ptr).u[M_A_Y]GRD_IND;

       m_b   = (*fields_ptr).u[M_B  ]GRD_IND;
       m_b_t = (*fields_ptr).u[M_B_T]GRD_IND;
       m_b_x = (*fields_ptr).u[M_B_X]GRD_IND;
       m_b_y = (*fields_ptr).u[M_B_Y]GRD_IND;

       m_c   = (*fields_ptr).u[M_C  ]GRD_IND;
       m_c_t = (*fields_ptr).u[M_C_T]GRD_IND;
       m_c_x = (*fields_ptr).u[M_C_X]GRD_IND;
       m_c_y = (*fields_ptr).u[M_C_Y]GRD_IND;

		m_a_up_x = g_up_xx * m_a_x + g_up_xy * m_a_y;
		m_a_up_y = g_up_xy * m_a_x + g_up_yy * m_a_y;
		
		m_b_up_x = g_up_xx * m_b_x + g_up_xy * m_b_y;
		m_b_up_y = g_up_xy * m_b_x + g_up_yy * m_b_y;

		m_c_up_x = g_up_xx * m_c_x + g_up_xy * m_c_y;
		m_c_up_y = g_up_xy * m_c_x + g_up_yy * m_c_y;

        norm_m_2 = (m_a*m_a + m_b*m_b + m_c*m_c);

		Eint = m_a_t*m_a_t + m_a_x*m_a_up_x + m_a_y*m_a_up_y 
		     + m_b_t*m_b_t + m_b_x*m_b_up_x + m_b_y*m_b_up_y
		     + m_c_t*m_c_t + m_c_x*m_c_up_x + m_c_y*m_c_up_y;
			  
		Eint = Eint -   (   (m_a*m_a_t + m_b*m_b_t + m_c*m_c_t)*(m_a*m_a_t + m_b*m_b_t + m_c*m_c_t)
					    +  (m_a*m_a_up_x + m_b*m_b_up_x + m_c*m_c_up_x)*(m_a*m_a_x + m_b*m_b_x + m_c*m_c_x)
		                +  (m_a*m_a_up_y + m_b*m_b_up_y + m_c*m_c_up_y)*(m_a*m_a_y + m_b*m_b_y + m_c*m_c_y)
		                )
		                /norm_m_2;

		E = E + Eint*Weight(grid_ptr,grid_ind1,grid_ind2)*detg_1_2;
     }
   }
 }
    E = E*one_dN1*one_dN2;
	

#ifdef MPI

{   int tag;
    int grid_to;
    int grid_from;
    int grid_ind;
    MPI_Status status;
    FLOAT EE[N_GRIDS];
    MPI_Request reqs[N_GRIDS],  reqr[N_GRIDS];

     	if(grid != 0){

		grid_to=0;
		tag=grid;
		MPI_Isend(&E,1, MPI_DOUBLE, grid_to, tag, global, &reqs[tag]);
		MPI_Wait(&reqs[tag], &status);
		E = 0.;
		}
	else {
		EE[0]=E;
     		for(grid_ind = 1; grid_ind < N_GRIDS; ++grid_ind){

     			grid_from=grid_ind;
     			tag=grid_ind;
     			MPI_Irecv(&EE[grid_ind],1, MPI_DOUBLE,  grid_from, tag, global, &reqr[tag]);
			MPI_Wait(&reqr[tag], &status);
                E = E + EE[grid_ind];
			}
		}
}
#endif // MPI

	return(E);
} // Function Energy




// ---------------------------------------------------------------------------------------------------------------------------------

FLOAT Winding_number(struct GRID_PAR *grid_ptr,
		struct field_array *fields_ptr, 
		struct field_array *derivs_ptr,
		struct FUNCTION_PAR *function_par) {

// Returns the winding number of a field, this is the integral of B_0 = -1/(4pi){m^x[d_x m^y d_y m^z - d_y m^y d_x m^z] - m^y[d_x m^x d_y m^z - d_y m^x d_x m^z] + m^z[d_x m^x d_y m^y - d_y m^x d_x m^y] }

    int ni_1 = (*grid_ptr).start_grid_1; 
    int nf_1 = (*grid_ptr).final_grid_1;  
    int ni_2 = (*grid_ptr).start_grid_2;  
    int nf_2 = (*grid_ptr).final_grid_2; 
    int n_gridpts_1 = nf_1 - ni_1;
    int n_gridpts_2 = nf_2 - ni_2;
    int grid = (*grid_ptr).grid;
    int coorx = 1;    
    int coory = 2;


    FLOAT xi = (*grid_ptr).initial_x;
    FLOAT xf = (*grid_ptr).final_x;
    FLOAT yi = (*grid_ptr).initial_y;
    FLOAT yf = (*grid_ptr).final_y;

    FLOAT x, y, m_t_nn, nn_x, nn_y, nn_z;

    FLOAT one_dN1 = (xf-xi)/(FLOAT)(nf_1-ni_1-1);
    FLOAT one_dN2 = (yf-yi)/(FLOAT)(nf_2-ni_2-1);

    FLOAT sigma = (*function_par).sigma; /* For dissipation */

    FLOAT mm = (*function_par).mm;
	FLOAT W = 0.0;
	
	FLOAT m_a=0.0;
    FLOAT m_a_x=0.0;
    FLOAT m_a_y=0.0;

    FLOAT m_b=0.0;
    FLOAT m_b_x=0.0;
    FLOAT m_b_y=0.0;

    FLOAT m_c=0.0;
    FLOAT m_c_x=0.0;
    FLOAT m_c_y=0.0;
	
//    FLOAT g_up_xx, g_up_yy, g_up_xy, g_up_yx, detg_1_2;


 {register int grid_ind1, grid_ind2;
 
 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){
   for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){

       m_a   = (*fields_ptr).u[M_A  ]GRD_IND;
       m_a_x = (*fields_ptr).u[M_A_X]GRD_IND;
       m_a_y = (*fields_ptr).u[M_A_Y]GRD_IND;

       m_b   = (*fields_ptr).u[M_B  ]GRD_IND;
       m_b_x = (*fields_ptr).u[M_B_X]GRD_IND;
       m_b_y = (*fields_ptr).u[M_B_Y]GRD_IND;

       m_c   = (*fields_ptr).u[M_C  ]GRD_IND;
       m_c_x = (*fields_ptr).u[M_C_X]GRD_IND;
       m_c_y = (*fields_ptr).u[M_C_Y]GRD_IND;


		W = W + Weight(grid_ptr,grid_ind1,grid_ind2)*( m_a * (m_b_x * m_c_y - m_b_y * m_c_x) 
		                                             - m_b * (m_a_x * m_c_y - m_a_y * m_c_x) 
		                                             + m_c * (m_a_x * m_b_y - m_a_y * m_b_x)
		                                             );
				  
     }
   }
 }
    W = -W*one_dN1*one_dN2/4./PI;
	

#ifdef MPI

{   int tag;
    int grid_to;
    int grid_from;
    int grid_ind;
    MPI_Status status;
    FLOAT EE[N_GRIDS];
    MPI_Request reqs[N_GRIDS],  reqr[N_GRIDS];

     	if(grid != 0){

		grid_to=0;
		tag=grid;
		MPI_Isend(&W,1, MPI_DOUBLE, grid_to, tag, global, &reqs[tag]);
		MPI_Wait(&reqs[tag], &status);
		W = 0.;
		}
	else {
		EE[0]=W;
     		for(grid_ind = 1; grid_ind < N_GRIDS; ++grid_ind){

     			grid_from=grid_ind;
     			tag=grid_ind;
     			MPI_Irecv(&EE[grid_ind],1, MPI_DOUBLE,  grid_from, tag, global, &reqr[tag]);
			MPI_Wait(&reqr[tag], &status);
                W = W + EE[grid_ind];
			}
		}
}
#endif // MPI

	return(W);
} // Function Winding_number


FLOAT Weight(struct GRID_PAR *grid_ptr,int i, int j) {
// computes the weights to make the integrals

    	int ni_1 = (*grid_ptr).start_grid_1; 
    	int nf_1 = (*grid_ptr).final_grid_1;  
    	int ni_2 = (*grid_ptr).start_grid_2;  
    	int nf_2 = (*grid_ptr).final_grid_2; 
	FLOAT W;

		if(((i==ni_1 || i==nf_1-1)  && (j!=ni_2 || j!=nf_2-1)) || ((i!=ni_1 || i!=nf_1-1)  && (j==ni_2 || j==nf_2-1))){W=0.5;}
			else{
				if((i==ni_1 || i==nf_1-1)  && (j==ni_2 || j==nf_2-1)){W=0.333333333333333;}
				else{W=1.;}
				}

return(W);
}
				
